Computational Vision
U. Minn. Psy 5036
Daniel Kersten
Lecture 4: Ideal Observer Analysis

Goals

Last time

The accuracy and reliability of human perceptual decisions are limited by two primary sources:

	1) inherent uncertainty in the stimulus information for a specific task
	2) limitations of the human observer.

Last time we focused on the notion of an ideal observer. The “ideal” is presumed to have knowledge of the generative process that produces the image data. Its ability to make reliable decisions is limited by 1).   

For the dot simulations of photon distributions, the variability is modeled as a Poisson distribution. We can think of a distribution as the limiting case of a histogram, where we have unlimited samples. The spread of the distribution is a measure of the variability in the data (dot count) for a fixed switch setting -- the variability in the generative process. This figure shows the distributions and histograms for high and low settings. These probability distributions of dot counts x, are said to be conditional on the switch setting, and we express this as: p (x|(Subscript[S, L])), p (x|(Subscript[S, H])).

nsamples =100;Show[Histogram[Table[RandomVariate[PoissonDistribution[i],nsamples],{i,{4,9}}],{1},"PDF",PlotRange->All,ChartStyle->{LightBlue,LightBrown}],Plot[Evaluate[Table[PDF[PoissonDistribution[μ],x],{μ,{4,7}}]],{x,0,15},PlotStyle->Thick,PlotRange->All,PlotLegends->{p(x|(Subscript[S, L])),p(x|(Subscript[S, H]))}]]

The ideal's job is to take input data and make the best guess regarding some state of the generating process--i.e. make an inference. In one demo, the image data was the display of dots, the task is to decide the switch state, high or low, and "best" is defined as minimizing the error rate. This is a called a yes/no discrimination task, where the observer is shown only one image in a given “trial”, but it could have been caused either by the high or low switch setting, and the observer has to decide which.  Intuitively, it seems like the best strategy is to count the number of dots, and use this number to decide the state. But how should one determine the criterion that has to be crossed before deciding the switch state is probably high vs. low?

In the case of a two-alternative forced-choice discrimination task, an observer is given two images in a given trial--one was caused by the high switch, and one by the low. The only source of information to make the decision is the image data--the dots (i.e. not which pattern was shown first, or whether shown on the left or right). The observer has to decide which of the two was caused by the "high" switch (or "low"). Whether the "signal" is defined as the high or low setting, is a matter of convention, and terminology. We'll treat the high switch setting as the "signal".

dotpattern[mean_]:=Table[RandomReal[{0,1},2],{RandomVariate[PoissonDistribution[mean]]}];display[mean_]:=Graphics[{PointSize[0.04],White,Point/@dotpattern[mean]},AspectRatio->1,Frame->False,FrameTicks->None,Background->GrayLevel[0.],PlotRange->{{-0.2,1.2},{-0.2,1.2}}];
Show[GraphicsRow[{display[9],display[4]},Spacings->Scaled[0.3]],Frame->False]

Appealing to intuition again suggests that the ideal strategy would be to count the dots in each, and decide “signal” for which ever one has the most dots. We’ll see that this indeed is the ideal strategy, but only under a particular generative model. For example, in the above case, you really don’t have to count the dots at all to make the best decision! This is because there is a “prior” condition, that puts the high switch display on the left 100% of the time. However, if we randomize which display comes first, you would have to use the data:

GraphicsRow[display[#]& /@RandomSample[{4,9},2]]

Today

In this lecture we complete our introduction to classical signal detection theory (SDT). SDT provides an important set of tools for measuring and modeling the sensitivity of human and neural perceptual decisions. (Later we'll generalize further to "statistical decision theory"--same acronym!) We will: 

	- Understand how to summarize ideal (and human performance) in the yes/no task in terms of hit and false alarm rates, and to relate these to a sensitivity measure called d’. To do this, we will introduce the (standard) Gaussian approximation, and apply it to variability in light levels.
	- Understand how to quantitatively estimate the “signal-to-noise” ratio “inside the head”, and from there compare human and ideal performance.
	- Calculate performance in other tasks. In particular, the two-alternative forced-choice task 
	- Measure your own absolute statistical efficiency in a 2AFC task

What we learn today will provide the basis for addressing the question: What does the eye see best?

Ideal inference model for dot density discrimination: yes/no task (from Lecture 3)

Our goal now is to show how the use of the dot count arises from a deeper principle:  The best guess should be based on choosing the most probable switch setting given a measurement. This is called choosing the setting that "maximizes the posterior probability".

Yes/no task

In a yes/no task, we designate one switch position as the "signal" and the other as "not the signal", or "noise". We'll have the observer say "yes" when the switch is set to high.

Posterior probability distribution

Given a measurement, what was the switch setting? We've seen that because of variability in the measurement, one can’t answer this question correctly every time just from the data. However, we can ask how the decision should be made so that an observer is correct as often as possible (equivalent to minimizing the probability of error). So we ask:

	 What decision strategy will minimize the probability of error?

The observer that achieves a pre-defined optimality criterion is our "ideal observer". We'll see other ideal criteria later.

Given a measurement,  assume that the ideal observer should decide “signal” (i.e. decide the hypothesis) of this hypothesis has the bigger probability. In other words, measure x, and then calculate:

If P((Subscript[S, L])| x) > P((Subscript[S, H])|x), choose (Subscript[S, L])
	otherwise, choose (Subscript[S, H]).

P((Subscript[S, L]) | x) is "the probability of the switch taking on the value of (Subscript[S, L]), conditional on a measurement x", and similarly for P((Subscript[S, H]) | x). 

This rule is called the "maximum a posteriori" or MAP rule. ...but these aren’t the conditional probabilities we calculated and graphed above! 

P((Subscript[S, L]) | x) and P((Subscript[S, H]) | x) are called posterior probabilities. Why "posterior"?  Even without any data measurements, knowledge of P((Subscript[S, L])) and P((Subscript[S, H])) can provide important information for optimal decisions, and these are called prior probabilities, i.e. the information the observer has prior to the arrival of the data measurements. Posterior probabilities can be thought of as prior probabilities that get updated by data.

Now it is often hard to directly specify the posterior probability, and it is easier to model the prior probabilities and the likelihoods, and then use Bayes rule to relate them to each other:

p(H|x)=((p(x|H)p(H))/(p(x)))∝p(x|S)p(S)

..ah, looks like now we will be able to plug in the conditional probabilites we calculated earlier. What we have done is to partition the generative components into two parts: a likelihood, p(x|H) and a prior, p(H).

Bayes' rule is a classic rule in probability theory, attributed to the Rev. Thomas Bayes (1702-1761). However, while Bayes introduced the idea of inverse probablistic inference, it was later, in the early part of the 18th century, that Laplace produced the above equation.  After you’ve gone over the basics of probability (see ProbabilityOverview.nb) you'll be able to easily prove it.

Optimal performance is determined by the posterior probability (and task), and that is why Bayes' rule is of central importance in determining optimal theories of inference. Let's consider the roles the priors and likelihoods play separately.

Prior only

Suppose the switch is set to (Subscript[S, L]) 4/5 of the time, and to (Subscript[S, H]) 1/5 of the time.

(*pH={1/2,1/2};*)
pH={4/5,1/5};

BarChart[pH,ChartStyle->{Red,Blue},ChartLegends->{"p(H=(Subscript[S, L]))","p(H=(Subscript[S, H]))"}]

But if the dot counter isn't working at all--there are no data. What is the observers's optimal strategy to guess whether it was the high or bright flash? 

If the receiver’s goal is to minimize the average probability of error, then it should always pick (Subscript[S, L]). Or to use the signal convention, always decide “signal not present”.

Like the maximum a posteriori rule, the rule that minimizes the error rate is:
	choose  (Subscript[S, H]) if p ((Subscript[S, H])) > p((Subscript[S, L]))
	choose (Subscript[S, L]) if p((Subscript[S, L])) > p( (Subscript[S, H]))

(Interesting enough, people often don't do this. Imagine you are the observer--wouldn’t it be really hard to say “no signal” every single time? With feedback, instead human observers and animals tend to match probabilities. This is a phenomenon which has been studied by cognitive psychologists and economists, cf. Green et al., 2010).

Following the above rule, (Subscript[S, H]) is presented p((Subscript[S, H])) fraction of the time--so p((Subscript[S, H])) = 1/5 is exactly the error rate, i.e. the fraction of the time that the observer gets the wrong answer using the minimum error rule.

Here is a simple simulation of the process, and the error rate after numtrials:

numtrials = 100000;
trials = RandomVariate[BernoulliDistribution[1/5],numtrials];(*Bernoulli sampling is just a coin flip--biased to favor 0 over 1 in our case*)
idealdecisions = Table[0,numtrials]; (*0 stands for always deciding "no signal"*)
Count[MapThread[Equal[#1,#2]&,{trials,idealdecisions}],False]/numtrials//N
Rationalize[%,.1]

Of course we could just calculate: Count[trials,0]/numtrials. But...

Exercise: what is the probability of error if the probabilities are matched? Here is a simulation of the process, but can you derive the answer?

probmatchdecisions = RandomVariate[BernoulliDistribution[1/5],numtrials]; (*randomly decide "0" 4/5th of the time*)

Count[MapThread[Equal[#1,#2]&,{trials,probmatchdecisions}],False]/numtrials //N(*count errors*)

Likelihood

We use Mathematica to define p(x|H) = poisson[x, a], where the hypothesis or switch setting H, determines the average number of samples, x:

poisson[x_,a_]:=PDF[PoissonDistribution[a],x]

Then we make two lists of probabilities corresponding to means of 4 and 9:

p1 = Table[N[poisson[x,4]], {x,0,20}];
p2= Table[N[poisson[x,9]], {x,0,20}];

BarChart[Transpose[{p1,p2}],PlotRange->{0.,.25},AxesLabel->{"x","p(x|H)"},ChartStyle->{Red,Blue},ChartLegends->{"p(x|H=(Subscript[S, L]))","p(x|H=(Subscript[S, H]))"},BarSpacing->-.55]

This is a barchart version of the plot we showed earlier, showing the frequencies of observing x dots under the two possible switch settings, H=(Subscript[S, L]) and H=(Subscript[S, H]).

Let's put the two likelihoods together into one 2x21 array:

poispxH = {p1,p2};
Dimensions[poispxH]

Combining the prior and likelihood: Maximum a posteriori (MAP) rule

By definition, the ideal “knows” the likelihoods: p(x|(Subscript[S, L])) and p(x|(Subscript[S, H])) and the priors, p((Subscript[S, L])) and p((Subscript[S, L])). 
Bayes' rule tells us how to combine them to obtain the a posteriori probability of the hypotheses conditional on the data x:

p(H|x)=((p(x|H)p(H))/(p(x))) =((p(x|H)p(H))/((∑_H)p(x|H)p(H)))

(The expression for the denominator p(x) follows from the basic laws of probability, see notebook / pdf.).

Let’s first asssume uniform priors:

pH={1/2,1/2};

(*pH={4/5,1/5};*)

For our example, we can calculate the posterior using the following, concise (if initially opaque), line of Mathematica code:

posterior = Transpose[Transpose[(poispxH pH)]/Plus@@(poispxH pH)];

Below is a plot of the posterior probabilities, p((Subscript[S, H])|x) and  p((Subscript[S, L])|x),  as function of the data x:

BarChart[Transpose[{posterior[[1]],posterior[[2]]}],PlotRange->{-.1,1.1},AxesLabel->{"x","p(H|x)"},ChartStyle->{Red,Blue},BarSpacing->-.55]

Note that the plot shows the probability of the switch being low or high, i.e. p(H=(Subscript[S, L])) or p(H=(Subscript[S, H])) given that x dots  were counted. 

(Why does the posterior probability plot look so different from the likelihood plot?)

Let's apply the maximum a posteriori rule (MAP)

For maximum a posteriori estimation, the rule is:

	choose (Subscript[S, H]) if p( (Subscript[S, H])|x) > p((Subscript[S, L])|x)
	choose (Subscript[S, L]) if p((Subscript[S, L])|x) > p( (Subscript[S, H])|x)

More generally, a function which chooses the value of argument (H = (Subscript[S, L]) or (Subscript[S, H])) that maximizes a function is sometimes written as:
 	(Underscript[argmax, H]) p(H|x)

Suppose 7 dots were counted, pick the most probable switch setting by inspecting the above graph. NOTE: the position of first bar corresponds to 0 dots.

Make a new plot with pH = {4/5, 1/5}. Now what is the most probable hypothesis?

Let's apply the Maximum likelihood rule

In maximum likelihood estimation, the rule is:

	choose (Subscript[S, H]) if p(x | (Subscript[S, H])) > p(x | (Subscript[S, L]))
	choose (Subscript[S, L]) if p(x | (Subscript[S, L])) > p(x | (Subscript[S, H]))

Following Bayes rule, if the prior probabilities are equal, then MAP is equivalent to maximum likelihood. 

p(x | (Subscript[S, H]))/p(x | (Subscript[S, L])) is called the likelihood ratio.

Suppose 7 dots were counted, what is the maximum likelihood decision? (See p(x|H) plots above)

Simple and classic OPTIMAL Rule for discrimination:  Compare the count to a fixed criterion.

Let's assume the prior probabilities are equal, so we base our decisions using the maximum likelihood rule. We will use Mathematica to do some symbolic manipulations to find the transition point in our data measurements (i.e. dot or photon count) where the likelihood goes from favoring a low to a high switch setting. This is where the likelihood ratio is 1, or the log of the ratio is zero.

loglikelihoodratio= Log[poisson[x1,b]/poisson[x1,d]]

PowerExpand[Log[(((b^x1) (E^(-b)))/(x1!))]-Log[(((d^x1) (E^(-d)))/(x1!))]]

There is a simple point to make here--the dot count, x1, is monotonically related to the log likelihood.  This means that our intuition that we can't do any better than simply counting dots was correct. But we need to decide on an appropriate criterion.

Solve[Log[(((b^x1) (E^(-b)))/(x1!))]==Log[(((d^x1) (E^(-d)))/(x1!))],x1]

If d and b are the average number of dots for the low and high switch settings, and the dot counts are distributed according to a Poisson distribution, then minimal error will result on average with the following rule:

 	Say "high" if 
		x>((b-d)/(Log[b]-Log[d]))
	
	Say "low" otherwise. 

The dot count, x, is said to be the decision variable. Because it is monotonically related to the likelihood ratio, it is an optimal decision variable. The transition value is called the  decision criterion: 

	decision criterion = ((b-d)/(Log[b]-Log[d]))
	
So now you have all you need to know to write a program that would make optimal decisions about high vs. low in a yes/no experiment. It would simply count dots, and base its decision on whether the count was greater or less than the decision criterion.

Summarizing performance using performance metrics

So far so good. But we would like to go a littler farther so that we can calculate theoretically what the average performance of an ideal observer would be.  

In a yes/no task, an ideal observer can make two kinds of mistakes. It can decide the signal was there when it wasn't (a false alarm also called a false positive), or it can decide it wasn't there when it was (a miss also called a false negative). These two determine the error rate. 

An ideal can also be correct in two ways. It can have a hit (true positive) or a correct rejection (true negative). We'll see below that you only need to measure two of these four statistics, because the hit + miss rate have to add up to 100% (if the signal is presented 100 times, you either say yes or no, so the sum of the hit and misses has to be 100).

We develop summary measures of discrimination performance in the next section.

The “standard” additive noise Gaussian model for signal detection

Motivation

Most inference modeling is done using Gaussian models of variability. One reason is theoretical convenience. A deeper theoretical reason rests on the Central Limit Theorem, which says that a sum of independently drawn random variables (from a fixed distribution, even non-Gaussian) looks more and more Gaussian the more elements that are in the sum. Empirically, many experiments on human signal detection have been well-fit by assuming Gaussian distributions of the underlying decision variable. However, as we will see later (when we measure statistics on natural images), the Gaussian assumption/approximation for some random variables is a bad approximation. It is always important to test this assumption. We'll first show that the Gaussian approximation provides a good approximation to the Poisson distribution.

Some terminology. We've adopted the convention of treating the high or brighter light as a "signal". Similarly, in the context of detection, we can think of the low switch settings as "noise". We will continue with this here, and use the terms "signal" and "noise". But remember that this is just a convention--the problem is symmetric, and we could be talking about whether a measurement is from hypothesis A vs. hypothesis B (“Not A”).

What does i.i.d. mean?

Gaussian approximation for signal and noise Poisson distributions

As the mean a gets large, the frequency of occurrance of a Poisson distributed random variable can be approximated very well by the Gaussian  distribution:

p(X=x)=((e^(-((((x-μ))^2)/(2 (σ^2)))))/((Sqrt[2 π]) σ)). The mean or expectation of X is:E(X)=μ,and the variance is:var(X) =(σ^2)

This approximation is useful to estimate probability values for large a. If a is large enough, the probability of negative values (which is meaningless for a Poisson distribution) is very small. For computational convenience and for later generality, we will usually use the Gaussian approximation.

Let's compare the forms of the Poisson and Gaussian distributions:

Manipulate[
ndist=NormalDistribution[mean,(Sqrt[mean])];
pdist=PoissonDistribution[mean];
p1=Table[{x,PDF[pdist,x]},{x,-5,50}];
g1=ListPlot[p1,PlotStyle->RGBColor[1,0,0]];
p2=Table[{x,PDF[ndist,x]},{x,-5,50}];
g2=ListPlot[p2,Joined->True,PlotRange->{{-5,50},{0,.3}}];
Show[{g1,g2},PlotRange->{{-5,50},{0,.3}}],{{mean,10},2,40,1}]

Try comparing the Poisson and Gaussian approximation when the mean is much smaller, like a=4.

Note: Discrete vs. continuous distributions

The Poisson distribution represents probabilities of a random variable taking on integer values. The distribution is said to be "discrete". In contrast, the Gaussian distribution is continuous. A continuous distribution is represented by a probability density  meaning that we use it to characterize the probability of a Gaussian random variable falling within a certain range of real numbers.

We can interpret this approximation in two ways. We can discretize the continuous Gaussian function (as above) to give us a set of  probabilities (over integers) that closely match those of the corresponding Poisson distribution, and make sure that the discrete sum is one (a fundamental requirement for a probability distribution). Alternatively, as the photon count gets high, we can treat light intensity as a continuous quantity (abandoning our quantized notion of light magnitude). In this latter case, we would treat the random variable X (light intensity) as being a continuous variable with a continuous probability distribution or "density". Then, because there is an infinite number of possible values over any finite range, the probability of X=x, for any particular value (x = ℼ, or x = 3.1, for example) is actually zero! 

To fix this, we treat p(X) as a density (as in mass density in physics), rather than a probability (as in mass).  Then we can put a non-zero number on the probability of X taking on a value x in some small region, dx as: p(x<X<x+dx) ~p(x)dx. More on this later.

The generative model for signal discrimination: additive gaussian noise

Let's approximate our photon inspired model with a view towards generalization. We will express the generative model as an "additive gaussian model". This is a standard form used for all kinds of detection tasks, for visual and auditory patterns, as well as non-perceptual decisions. We can model the shift of the peak of the distribution as an additive offset to the mean of a Gaussian. Then we have:

H = (Subscript[S, H]): x=  b + noise;
	H = (Subscript[S, L]): x = d + noise;
	
where noise is a Gaussian distributed random variable with mean, μ = 0, and standard deviation σ. 
For the photon counting case, b=highmean, and d=lowmean. ("b for bright" and "d for dim"). Note that the standard deviations of the high and low distributions would, for a Poisson distribution, be different (variance = mean for Poisson). We will assume that for a typical discrimination task, the distributions are quite close together, so the standard deviations are almost equal. The assumption of Gaussian distributions with equal variance is common, because it simplifies calculations, but more importantly because in many practical cases of  discrimination, the approximation is pretty good.

Here is a plot of the theoretically predicted histograms for a signal (high) mean of 15, a noise (low) mean of 10, and a standard deviation of 4 for each:

gauss[x_,mean_,std_]:=((E^(-((((x-mean))^2)/(2 (std^2)))))/(std (Sqrt[2 π]))); (*Define our own gaussian distribution*)
b=15;d=10;sigma=4; max=gauss[0,0,sigma];
Manipulate[
Plot[{gauss[x,d,sigma2],gauss[x,b1,sigma2]},{x,-5,30},AxesLabel->{"x","p"},Filling->Axis,PlotRange->{0,max +0.25},Epilog->{Text["(Subscript[μ, S])=b",{b1,0.11}],Text["σ",{b1+sigma2*1.4,(Exp[-.5]/(Sqrt[2.0*Pi]*sigma2))}],Line[{{b1,(Exp[-.5]/(Sqrt[2.0*Pi]*sigma2))},{b1+sigma2,(Exp[-.5]/(Sqrt[2.0*Pi]*sigma2))}}],Text["(Subscript[μ, N])=d",{d,0.11}]}],{{b1,b, "μ=b"},d,30},{{sigma2,sigma, "σ"},1,6}]

The signal-to-noise ratio: d', a summary statistic for ideal performance

Whether a measurement x comes from the signal or the noise is less ambiguous if the difference in means, b-d, is big. But it is also less ambiguous if the standard deviation σ is smaller. By using Gaussian distributions (with equal variances), we can characterize the ideal's signal-to-noise ratio with one number, the "signal-to-noise ratio," defined as  d’:

		d'=((b-d)/σ)
		
where b and d are the high and low means, respectively.  This makes intuitive sense. Discrimination should get easier as the difference between the means increases (the "signal" is the difference) or as the spread given by the standard deviation of the additive noise decreases--hence the term signal-to-noise ratio. Note that d' can also decrease by increasing σ.

But what does this mean in terms of performance?

The inference model for signal discrimination given additive gaussian noise

How does the ideal observer make a decision as to whether the low or high light was flashed? Earlier we derived a criterion starting from the assumption that we wanted to  maximize the posterior probability (p(H|x) over H:

1) 	(Underscript[argmax, H]) p(H|x)

Which means "find that value of H (e.g. H = “switch high” vs. H= “switch low”) which makes p(H|x) the biggest".  From here, we showed that if the prior probabilities over the hypotheses were the same (p(H=(Subscript[S, H])) = p(H=(Subscript[S, L])), this was equivalent to  maximizing the likelihood:

2) 	(Underscript[argmax, H]) p(x|H)
Because we are considering only two hypotheses, we could reformulate the decision strategy to testing whether the ratio, ((p(x|(Subscript[S, H])))/(p(x|(Subscript[S, L]))))>1 ? This in turn, is equivalent to testing: log[((p(x|(Subscript[S, H])))/(p(x|(Subscript[S, L]))))]>0? 

Earlier we applied the likelihod ratio rule to our dot density version of light intensity discrimination. We showed the decision could be based on whether the photon count was bigger than a particular criterion (call it (Subscript[X, T])) determined by two light level means (b and d).

Graphically, the ideal that minimizes its error rate makes its decision by deciding "high" if the measurement x is right of the cross-over point on the above plot (i.e. x where  ((p(x|(Subscript[S, H])))/(p(x|(Subscript[S, L])))) = 1 )

This minimizes the probability of error, but how is error related to the decision criterion and the distributions?

Let's consider the more general case, where the criterion isn't necessarily at the cross-over point.

Hit, false alarm (positive), miss, and correct rejection rates.

It is easy to imagine how one might  experimentally measure the signal-to-noise ratio for light discrimination for the ideal observer—we just collect histograms under the two conditions (H=(Subscript[S, L]) and H = (Subscript[S, H])), approximate them by Gaussian distributions (giving us two conditional probability distributions), and assuming the standard deviations are close, use these Gaussian fits to estimate d'. We simulated doing this kind of thing in the last couple of lectures.

But how could we possibily measure the signal-to-noise ratio, or d' of a human observer?
 
Human decisions are based on some hidden, and probably quite complex neural mechanism in the brain (cf. Yang & Shadlen, 2007). It seems like we'd need to have access to a neural response that behaves like the ideal's decision variable, but is consistent with human performance (which is usually sub-ideal). This is an interesting scientific problem, but let's see if we can put a number on human d', without "going inside the box"--just based on behavioral performance.

To answer this question, let's take a look at an alternative way of estimating d' for the ideal observer. The ideal observer (or receiver) for light intensity discrimination has two ways of being right and two ways of being wrong:

Being correct:
	it can score a
		 hit (e.g. says "high" when the switch was set to high)
	or a
		correct rejection (e.g. says "low" when the switch was set to low)
Being incorrect:
	it can suffer a 
		false alarm (also called "false positive") (e.g. says "high" when the switch was set to low)
	or a 
		miss (or "false negative") (e.g. says "low" when the switch was set to high)
(Statisticians talk about a similar distinction in terms of Type I (false positive) and Type II (false negative) errors).

Rates

Average performance in a yes/no task is completely characterized by calculating the proportions of two of the four. Hit and false alarm rates can be treated as estimates of conditional probability distributions, p(response | switch setting, H). For example, 

	hit rate =  ((#times observer says "high when switch was set to high")/(#times switch was set to high)) 
	 ~ p(decide high | switch set on high)
	 
	false alarm (positive) rate =   ((#times observer says "high"when switch was set to low)/(#times switch was set to low)) 
 	~ p(decide high | switch is set to low)

	Sometimes, we talk about the average error rate. Since there are two ways of being wrong: Deciding "high" when H = (Subscript[S, L]), and deciding "low", when H = (Subscript[S, H]). The total error rate is the (weighted) average of the miss and false alarm rates. The error rate is determined by the mean values for the high and low settings. As b increases, the separation between the probability distributions increases, and the overlap decreases, so the error rate decreases. So intuitively, there should be some relationship between d' and error and/or success rates.

We only need measures of hit and false alarm rates because the they are related to the correction rejection and miss rates.  Show that:

The corresponding correct rejection and miss rates are: 

	p(correct rejection)  = 1 - p(false alarm), 
and 
	p(miss) = 1 - p(hit),  respectively.

How should one compute weighted average for error rate?

Graphical view of  the hit, false alarm rate, ...

For a probability density (continuous distribution) function (i.e. a "PDF"), say p(x), the probability of a measurement X falling within a certain range is given by the area under the density over that range:

	P(x1 <X < x2) = (∫_x1^x2)p(x)dx
	P(X>XT) = (∫_XT^∞)p(x)dx

The criterion, and thus the hit and false alarm rates could be determined by the relative costs or benefits (loss or gain) one assigns to a particular choice of hit and false alarm rates. Suppose the criterion is XT, in general not at the cross-over point of the likelihoods (which would only be optimal for constant prior probabilities and the goal of minimizing average error). Then the hit rate (PH) is determined by the area under the  (signal or high) curve to the right of XT. The false alarm rate (PFA) is given by the area under the ("noise" or low) curve to the right of XT.

b=15; d=10; sigma = 4;max = 1.1;
ndistd=NormalDistribution[d,sigma];
ndistb=NormalDistribution[b,sigma];
max=PDF[ndistb,b];
Manipulate[
g1=Plot[{PDF[ndistd,x],PDF[ndistb,x],(UnitStep[x-c]*Max[PDF[ndistb,x],PDF[ndistd,x]])},{x,-5,30},AxesLabel->{"x","p"},Filling->Axis,PlotRange->{0,max +0.025},Epilog->{Text[Style["μ=b",7],{b,0.11}],Text[Style["μ=d",7],{d,0.11}]}];
g2=Plot[{1-CDF[ndistd,x],1-CDF[ndistb,x],(UnitStep[x-c]*(1-CDF[ndistb,x]))},{x,-5,30},AxesLabel->{"x","Hit & FA"}];
GraphicsGrid[{{g1,g2}}],{{c,b,"criterion"},0,30}]

Criterion shifts affect hit and false alarm rates, but don't affect basic sensitivity

Consider again our decision rule:
	if x > (Subscript[X, T]) guess "high switch caused the intensity measured"
	if x<=(Subscript[X, T]) guess "low switch caused the intensity measured"
In general, where the criterion gets placed depends on the decision goal.  One could have other  goals (than minimizing error) that would determine where to put the criterion level.  Put yourself in the place of an ideal (not a MAP observer) with the following constraints:
	If you were slapped on the wrist every time you said "high", you might never say high--you would never get any hits. This in effect pushes the criterion far to the right.
	If you liked chocolates as much as I do, and received a sweet every time you said high, you might always say high, even if you thought the signal  was not presented-- after all, why not be optimistic? You would have many false alarms. This pushes the criterion far  to the left.
	
So the goal doesn't have to be determined by maximizing the proportion of correct responses (minimizing error), it can be determined by other criteria, which in turn modify the decision rule. These other factors can be incorporated into a pay-off matrix (see Green and Swets, 1974). 

In that we haven't changed the means or standard deviations, d' hasn't changed. Hit and correct rejection rates trade-off against each other. You can get more hits but at the expense of making more false alarm mistakes. (Recall correction rejection rate = 1 - hit rate). It seems like there should be some way to calculate d' from the hit and false alarm rates. We'll see how to do that shortly. Later we will look at formalizing and generalizing the notions of costs and benefits as statistical decision theory. For example, the cost to an error in an estimate of illumination can be low as compared to a cost in the error of face identification.

Summary of modified log likelihood rule:

Thus we can see that one could derive a simple modification of the log likelihood rule:
Rather than testing: 
	log[((p(x|(Subscript[S, H])))/(p(x|(Subscript[S, L]))))] >  0?, 
instead we decide using:
	 log[((p(x|(Subscript[S, H])))/(p(x|(Subscript[S, L]))))]> k?, (equivalent to ((p(x|(Subscript[S, H])))/(p(x|(Subscript[S, L]))))> (e^k) )
where k is a function of the costs and benefits. Optimal decisions are based on the value of likelihood ratio. This ratio (or any monotonic function of it) is called the decision variable. The photon (or dot) count is a decision variable. More generally, a decision variable may or may not lead to optimal performance. For example, you may base your decision on a dot count, but you might have missed some.

Relationship between signal-to-noise ratio (d') and hit and false alarm rates.

We are now ready to see how to estimate the signal-to-noise ratio "inside an observer's head" using only performance measures of hit and false alarm rates. We need this to compare human d’ with an ideal d’. 

We noted that the signal-to-noise ratio d' can be estimated from the means and standard deviation. But if we don't have access to those numbers (as happens in psychophysics), we need a method to go from performance (e.g. hit and false alarm rates) back to a function of a hypothetical mean and standard deviation.

It turns out that with a bit of mathematics, one can show that d' can be obtained from the hit and false alarm rates using the following formula:

d'=z(PFA)-z(PH)

where z(p) is the z-score given a probability p. The figure below illustrates the relationship between probability P and z using the standard normal density (i.e. gaussian with zero mean and a standard deviation of 1):

P=(∫_z^∞)gauss[x,0,1]dx. And z(p) is the inverse.

There is no simple formula for z, but there are good closed-form approximations. Mathematica doesn't give the direct formula for the z-score, but it does supply the inverse of the erf[] function, which comes from the engineering (rather than statistics) tradition:

erf(z)=(2/(Sqrt[π]))(∫_0^z)(e^(-(t^2)))dt

The z-score function is related to the InverseErf function by:

z[p_]:=Sqrt[2] InverseErf[1-2 p];

You will use the z-score function in Assignment #1.

Statistical efficiency of human vision: Humans vs. ideals

The basic idea: human observers are sub-ideal, but may be "ideal-like"

The idea is to model the human signal discrimination  as being "ideal-like" in assuming that human decisions respect an implicit generative model:

	H = (Subscript[S, H]): x= (b^h) + noise';
	H = (Subscript[S, L]): x = (d^h)+ noise';
	
where (b^h), (d^h), and noise' are the equivalent states of the world (corresponding to the two effective means of the human observer) that could give rise to the human's d', as measured from hit and false alarm rates:

(Subscript[d', human])=z((Subscript[PFA, human]))-z((Subscript[PH, human]))

It is as if the human visual system does Bayesian inference (is ideal), but has the wrong generative model--i.e. a different state of the world. The d' for human, is determined by hit and false alarm rates, or equivalently by: 

	d'for human=(((b^h)-(d^h))/(σ^h)), 
Note that there is indeterminancy in these "implicit" variables, (b^h), (d^h), and (σ^h) (the standard deviation of noise') -- there is an infinite family of combinations of (b^h), (d^h), and (σ^h)  which give the same d'.
How good is the model? One way of testing it is to plot hit and false alarm rates for human decisions and compare them to this "sub-ideal" that has additive gaussian noise with equal variances. Surprisingly often, the model fits are quite good. But first, lets see how we can make an absolute comparison of performance.

Comparing ideal and human performance

For the light discrimination problem, the physics of the experiment determines the generative model, i.e. the mean levels  d, b and the standard deviation. We have seen that the ideal's performance is characterized by one number called the sensitivity d'. Now that we understand the limitations on the performance of an ideal observer, let’s  how to compare human performance to the ideal. Even if the ideal is making near perfect discriminations, the human observer may not be doing so well because of other sources of uncertainty. For example, the ideal may be contending with the following situation:

We can't "see" or directly measure the distributions that the human observer is using to make the decision, but we can suppose that it is based on distributions that are in effect much closer together:

Or they could be noisier--i.e. bigger standard deviation than the ideal is coping with. 
Although we can not measure the human's d' by measuring the separation between these two distributions and their standard deviations (they are not directly measurable), we've seen that the human (or ideal) observer's performance can be determined by the hit and false alarm rate: d'for human =z(false alarm rate for human)-z( hit rate for human)

Statistical efficiency

Given the means to compute d' for the ideal and for the human observer in the same task, we can compare them. Usually we calculate the ideal's d' from the signal-to-noise ratio, and the human's from the hit and false alarm rate in a yes/no task or from the proportion correct in a 2AFC task (as in the experiment below). (The ideal's d' could be calculated from its hit and false alarm rates but this usually isn't as convenient--but it might serve as a good way to double-check that you are doing the right calculations.)

With these two d's in hand we can compare the performances of the two observers. One way is in terms of statistical efficiency.  Efficiency can be defined as the ratio (usually squared) of the two signal-to-noise ratios above:

	Statistical efficiency = (((d' for human/d' for ideal))^2). *

In our case of photon or dot density discrimination, one can interpret this ratio as the number of samples (e.g. photons) required by the ideal divided by the number of samples required by the human, when they are performing equivalently (e.g. same hit and false alarm rates). 

*(It is the reciprocal of this if the d' represents the physical signal to noise ratios at threshold. See Kersten and Mamassian, 2008)

Historical note -- Quantum Efficiency  accounting for the missing  information

In 1962, Horace Barlow reported results on the measurements quantum efficiency for light discrimination (rather than detection) under low light (scotopic) conditions similar to those of Hecht et al., and came up with a figure for QE of about 10%.  That is, the human observer behaved like an ideal observer who was only receiving one out of every ten photons. Where was the photon loss coming from? Like we saw for Hecht et al., Barlow traced the losses to reflection, scatter and absorption by the optic media,  and losses due to photons falling in the spaces between the rods, and an imperfect isomerization efficiency. Recall that a figure of 10% is close to what one would predict from Hecht et al.'s experiment.
	Barlow later went one step beyond Hecht et al.. He concluded (Barlow, 1977) that there was still a residual inefficiency even after taking into account all the above causes, which he calculated as accounting for only 80% of the photon loss. He was left with about 50% of human discrimination efficiency due to limitations in the brain's ability to "count" point events. That is, for example, if 100 photons are incident on at the cornea of the eye, about 20 of these are reliably transduced and this information is sent to the brain. But he argued, the brain deals with this average of 20 photons with 50% efficiency--that is, the ideal's "brain" could discriminate just as well with only an average of 10 photons. Barlow made this latter conclusion by a clever argument involving  a psychophysical experiment in which he had observers discriminate differences in dot density (rather than photon density) on a CRT screen. The idea was that although the presence of a photon at the retina does not necessarily make it to the brain, a dot will.

Psychophysical tasks & techniques

Testing our assumptions: The Receiver Operating Characteristic (ROC)

Although we can't directly measure the internal distributions of a human observer's decision variable, we've seen that we can measure hit and false alarm rates, and thus d'. But one can do more, and actually test to see if an observer's decisions are consistent with Gaussian distributions with equal variance. If the criterion is varied, we can obtain a set of n data points: 

{(hit rate 1, false alarm rate 1), (hit rate 2, false alarm rate 2), ..., (hit rate n, false alarm rate n)} 

all from one experimental condition (i.e. from one signal-to-noise ratio, call it (Subscript[d, ideal])'). This is because as the hit rate varies, so does the false alarm rate (see the above figures showing how hit and false alarm rates relate to area under the signal and noise distributions.). One could compute the d' for each pair and they should all be equal for the ideal observer. Of course, we would have to make a large number of measurements for each one--but on average, they should all be equal. 

	To get meaningful and equal d's for each pair of hit and false alarm rates assumes that the underlying relative separation of the signal and noise distributions remain unchanged and that the distributions are Gaussian, with equal standard deviation.  We might know this is true (or true to a good approximation) for the ideal, but we have no guarantee for the human observer. Is there a way to check? Suppose the signal and noise distributions look like:

If we plot the hit rate vs. false alarm rate data on a graph as the criterion xc varies, we get something that looks like:

b=15;d=10;sigma=4;
ndistd=NormalDistribution[d,sigma];ndistb=NormalDistribution[b,sigma];
max=PDF[ndistb,b];

Manipulate[
g1=Plot[{PDF[ndistd,x],PDF[ndistb,x],(UnitStep[x-c]*Max[PDF[ndistb,x],PDF[ndistd,x]])},{x,-5,30},AxesLabel->{"x","p"},Filling->Axis,PlotRange->{0,max +0.025},Epilog->{Text["μ=b",{b,0.11`}],Text["μ=d",{d,0.11`}]}];
g2=Plot[{1-CDF[ndistd,x],1-CDF[ndistb,x],(UnitStep[x-c]*(1-CDF[ndistb,x]))},{x,-5,30}];
g3=ParametricPlot[{{1-CDF[ndistd,x],1-CDF[ndistb,x]}},{x,-100,100},FrameLabel->{{"Hits",""},{"False Alarms","ROC curve"}},PlotRange->{{0,1},{0,1}},Frame->True,AspectRatio->1,Epilog->{Point[{1-CDF[ndistd,c],1-CDF[ndistb,c]}]}];
GraphicsGrid[{{g1,g2,g3}},ImageSize->Large],{{c,b,"criterion"},0,30}]

One can show that:

	 the area under the ROC curve is equal to the proportion correct in a two-alternative forced-choice experiment (see Green and Swets). 

From a purely empirical standpoint it can be useful to measure sensitivity by the area under the ROC curve. This provides a single summary number, even if the standard definition of d' is inappropriate, for example because the variances are not equal, or the distributions are not gaussian.

We return to our basic question: is there a way to spot whether our gaussian equal-variance assumptions are correct for human observers? If we take the same data and plot it in terms of Z-scores we transform the ROC curve to a straight line:

b=15;d=10;sigma=4;
ndistd=NormalDistribution[d,sigma];ndistb=NormalDistribution[b,sigma];
max=PDF[ndistb,b];
zscore[p_]:=(Sqrt[2]) InverseErf[1-2 p];

Manipulate[

g3=ParametricPlot[{{1-CDF[ndistd,x],1-CDF[ndistb,x]}},{x,-100,100},FrameLabel->{{"Hits",""},{"False Alarms","ROC curve"}},PlotRange->{{0,1},{0,1}},Frame->True,AspectRatio->1,Epilog->{Point[{1-CDF[ndistd,c],1-CDF[ndistb,c]}]}];

g4=ParametricPlot[{{zscore[1-CDF[ndistd,x]],zscore[1-CDF[ndistb,x]]}},{x,-100,100},FrameLabel->{{"z(Hits)",""},{"z(False Alarms)","Z score ROC curve"}},PlotRange->{{-10,10},{-10,10}},Frame->True,AspectRatio->1,Epilog->{Point[{zscore[1-CDF[ndistd,c]],zscore[1-CDF[ndistb,c]]}]}];

GraphicsGrid[{{g3,g4}}],{{c,b,"criterion"},0,30}]

In fact, if the underlying distributions are Gaussian, the data should lie on a straight-line. If they both have equal variance, the slope of the line should be equal to one. This is because:

And if we solve for the criterion Xc, we obtain:

(I've switched notation here, where b = μs, and d = μn). The main point of this plot is to see if the data tend to fall on a straight line with slope of one. If a straight line, this would support the Gaussian assumption. A slope = 1 supports the assumption of equal variance Gaussian distributions.

	In practice, there are several ways of obtaining an ROC curve in human psychophysical experiments. One can vary the criterion that an observer adopts by varying the proportion of times the signal is presented. As observers get used to the signal being presented, for example, 80% of the time, they become biased to assume the signal is present. One needs to block trials in groups of, say 400 trials per block, where the signal and noise priors are fixed for a given block.
	
	One can also use a rating scale method in which the observer is asked to say how confident she/he was (e.g. 5 definitely, 4 quite probable, 3 don't know for sure, 2, unlikely, 1 definitely not). Then we can bin the proportion of "5's" when the signal vs. noise was present to calculate hit and false alarm rates for that rating, do the same for the "4's", and so forth. The assumption is that an observer can maintain not just one stable criterion, but four---the observer has in effect divided up the decision variable (x) domain  into 5 regions. An advantage of the rating scale method is efficiency--relatively few trials are required to get an ROC curve. Further, in some experiments, ratings seem psychologically natural to make. But if there is any "noise" in the decision criterion itself, e.g. due to memory drift, or whatever, this will act to decrease the estimate of d' in both yes/no and rating methods.
	
	Usually rather than manipulating the criterion, we would rather do the experiment in such a way that it does not change. Is there a way to reduce the problem of a fluctuating criterion?

The 2AFC (two-alternative forced-choice) method

Relating performance (proportion correct) to signal-to-noise ratio, d'.

In psychophysics, the most common way to minimize the problem of a varying criterion is to use a two-alternative forced-choice procedure (2AFC). In a 2AFC task the observer is presented on each trial a pair of stimuli. One stimulus has the signal (e.g. high flash), and the other the noise (e.g. low flash). The order, however, is randomized. So if they are presented temporally, the signal or the noise might come first, but the observer doesn't know which from trial to trial. In the spatial version, the signal could be on the left of the computer screen with the noise on the right, or vice versa. One can show that for 2AFC:

d'= - (Sqrt[2 ])z(proportion correct)

Exercise: Prove d'= - (Sqrt[2 ])z(proportion correct)

If you want to prove this for yourself, here are a couple of hints--actually, a lot of hints. Let us imagine we are giving the light discrimination task to the ideal observer. We have two possibilities for signal presentation: Either the signal is on the left and the noise on the right, or the signal is on the right and the noise on the left. There are two ways of being right. The observer could say "on the left" when the signal is on the left, or "on the right" when the signal is on the right. For example, for the light detection experiment, a reasonable guess is that all the ideal observer would have to do is to count the number of photons on the left side of the screen and count the number on the right too. If the number on the left is bigger than the number on the right, the observer should say that the signal was on the left. Thus, a 2AFC decision variable would be the difference between the left and right decision variables, where each of these  is what we calculated for the yes/no experiment.

r = (Subscript[r, L])-(Subscript[r, R])

For example as you will see in Assignment 1, (Subscript[r, L]) and (Subscript[r, R]) for the SKE observer would be the dot products of the signal pattern template with observation image vectors on the left and right sides.
So, the probabililty of being correct is: 
pc = p(r>0|signal on left) p(signal on left) + p(r<0| signal on right)p(signal on right)

What is the probability distribution of r? Well, from probability rules (see: ProbabilityOverview.nb), 
	average(r) = μ2 - μ1 = (Subscript[μ, s]) - (Subscript[μ, n])
	var(r) = var((Subscript[r, L]))+var((Subscript[r, R])), so (Subscript[σ, r])  = (Sqrt[2])(Subscript[σ, (Subscript[r, L])])=(Sqrt[2])(Subscript[σ, (Subscript[r, R])])
(Because the mean of the sum of two independent random variables is the sum of their means and that the variance of the sum is the sum of the variances.)
If the signal is equally likely to appear on the left or the right,the probability of being correct is the area under the curve to the right of zero of the distribution of r:

(Note in the abov e figure: r = (Subscript[r, L])-(Subscript[r, R])=(Subscript[x, L])-(Subscript[X, R]),and  μ2-μ1=(Subscript[μ, s])-(Subscript[μ, n])=b-d)

Calculating the Pattern Ideal's d' for a two-alternative forced-choice experiment from a z-score of the proportion correct. (see Homework Assignment #1)

For our 2AFC experiment, the observer gets two images to compare. One has the signal plus noise, and the other just noise. But the observer doesn't know which one is which. An ideal strategy is to compute the cross-correlation decision variable for each image (i.e. the dot product between each image vector and an exact template of the signal one is looking for), and pick the image which gives the larger cross-correlation. This strategy will result in a single measureable number, the proportion correct, Pc.

d' for a 2AFC task is given by the formula:

Where as before, the Z-score can be calculated from the inverse of a standard mathematical function called Erf[] to get Z from a measured P.

z[p_] := Sqrt[2] InverseErf[1 - 2 p];

where Z(*) is the z-score for Pc , the proportion correct. And then,

dprime[x_] := N[-Sqrt[2] z[x]]

2AFC demo for dot density discrimination

Next time

Probability overview: Look at the notebook / pdf

From pixels to patterns: What does the eye see best?

Computing the ideal observer for patterns

Comparing psychophysical performance for pattern detection with properties of visual neurons in the brain

Assignment #1.

Notes

Cumulative distributions

The cumulative distribution gives the probability that the detector signals x<k photons. It is obtained by adding up the probabilities for all values less than k. For the cumulative density function, we integrate over all values less than k.   Here is the cumulative distribution for the discrete Poisson distribution with a mean of 20:

Plot[CDF[pdist,x],{x,10,40}]

What is the probability of detecting 50 or less photons when the mean is 20? It is virtually certain-- as you can see from the graph, the probability is almost 1. Here is the plot of the continuous normal distribution with a mean of 20, and a standard deviation of Sqrt[20]:

Plot[CDF[ndist,x],{x,10,40}]

References

Applebaum, D. (1996). Probability and Information. Cambridge, UK: Cambridge University Press.
Barlow, H. B. (1962). A method of determining the overall quantum efficiency of visual discriminations. Journal of Physiology (London), 160, 155-168.
Barlow, H. B. (1977). Retinal and central factors in human vision limited by noise. In B. H. B., & F. P. (Ed.), Photoreception in Vertebrates Academic Press.
Barlow, H. B., & Levick, W. R. (1969). Three factors limiting the reliable detection of light by retinal ganglion cells of the cat. J Physiol, 200(1), 1-24.
Duda, R. O., & Hart, P. E. (1973). Pattern classification and scene  analysis . New York.: John Wiley & Sons.
Jason M. Gold, Craig Abbey, Bosco S. Tjan, and Daniel Kersten (2009) Ideal Observers and Efficiency: Commemorating 50 Years of Tanner and Birdsall: Introduction. JOSA A, Vol. 26, Issue 11, pp. IO1-IO2        doi:10.1364/JOSAA.26.000IO1
Geisler, W. (1989). Sequential Ideal-Observer analysis of visual discriminations. Psychological Review, 96(2), 267-314.
Green, D. M., & Swets, J. A. (1974). Signal Detection Theory and Psychophysics . Huntington, New York: Robert E. Krieger Publishing Company.
Yang, T., & Shadlen, M. N. (2007). Probabilistic reasoning by neurons. Nature, 447(7148), 1075–1080. http://doi.org/10.1038/nature05852
Treutwein, B. (1993). Adaptive psychophysical procedures: A review. Vision Research
Van Trees, H. L. (1968). Detection, Estimation and Modulation Theory . New York: John Wiley and Sons.
Watson, A. B., & Pelli, D. G. (1983). QUEST: A Bayesian adaptive psychometric method. 33, 113-120.
Watson, Andrew B.  & Fitzhugh, A., (1990) The method of constant stimuli is inefficient  Perception & Psychophysics 47(1), 87-91.